%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------
close all
%----------------------------------------------------------------
% 1. Endogenous variables (in total 97 variables)
%----------------------------------------------------------------
var 

% Utility (5)
uc Euc m Em mh

% Allocation (7)
y c cc cl i k l

% Prices (9)
rk w pi yn1 ph pe pih pie pihnum

% Multipliers (2)
lambda qk

% Productivity (2)
a dn 

% volatility (1)
sigmaa

% Energy sector (4)
eta ce epe etalag

% Growth (9)
dc dcl dcc dce dcr dz Edz dy di

% Auxiliary (9)
Phik Phik1 Fk Fl Jk Jl Phip Phip1 Phip2

% asset prices (49)
rf d dagg dd ddn ddr ddagg ddaggn ddaggr pd pdagg rd rdr rdagg rdaggn rdaggr rm rmex yn2 yn3 yn4 yn5 yn6 yn7 yn8 yn9 yn10 yn11 yn12 yn13 yn14 yn15 
yn16 yn17 yn18 yn19 yn20 slope fen cf dr cfn drn drrf dry rmagg rmexagg Erd Erdn 
% mf rf_for
;

%----------------------------------------------------------------
% 2. Exogenous variables
%----------------------------------------------------------------
 
varexo 
en em ee eeta ed esigma 
% efx
;

%----------------------------------------------------------------
% 3. Parameters
%----------------------------------------------------------------

parameters 

% Preferences
delta psi gamma xil omega lbar l_ss mu alphac phic

% Technology
nu alpha deltak phip xik rhon phin a1k a2k philev rhopd pd_ss phia

% productivity
rhoa a0

% volatility
rhosigma sigmaa0 sigmasigma sigmatarget rhosigma_a corr_vol

% dividend
sigmad

% Energy
sigmae sigmaeta rhoe rhoeta1 rhoeta2 deltaa ecap

% policy
rhor rhopi rhoy sigmam pi_ss yn1_ss 

% foreign sdf
% rhofx xifx

% steady state parameters
y_ss cc_ss ce_ss c_ss k_ss i_ss m_ss r_ss rf_ss 
Fl_ss Fk_ss Jl_ss Jk_ss w_ss qk_ss pe_ss Phik_ss Phik1_ss ph_ss;

%----------------------------------------------------------------
% 4. Calibration
%----------------------------------------------------------------
% Preferences
delta  = 0.994;
psi    = 2;
gamma  = 10;
lbar = 1;
alphac = 0.9;
phic = 5; % 5, 0.8 for Table 12 (3)
xil = 0.5;

%% Technology
nu = 6; % 6, 3 for Table 12 (2) and Table 13 (5)(6)(7)
alpha = 0.33;
deltak = 0.02;
phip = 80; % 80,30 For Table 12 (1) and Table 13 (4)(6)(7)
xik = 4.8;
rhon = 0.95; % 0.95, 0.9 for Table 12 (4)

% Productivity
rhoa = 0.983; % 0.983
a0 = -0.4471;
mu = 0.0046;
phia = 1; % 1, 2/3 for Table 13 (7)
phin = 0.1*sqrt(1-0.983^2)/(sqrt(1-rhoa^2))*sqrt(1-rhon^2)/sqrt(1-0.95^2); 
sigmaa0 = 0.018/sqrt(1-0.983^2)*sqrt(1-rhoa^2);

% Volatility
rhosigma = 0.997;
sigmatarget = 0.005;
corr_vol = 1;
rhosigma_a = corr_vol*sigmatarget/sigmaa0;
sigmasigma = sqrt(max(sigmatarget^2-rhosigma_a^2*sigmaa0^2,0));

% Dividend
sigmad = 0.08;

% Energy
rhoe = 0.6;
rhoeta1 = 0.99;
rhoeta2 = 0;
sigmae = 0.028*sqrt(1-rhoe^2); % 0.028*sqrt(1-rhoe^2) for benchmark, 0.2*sqrt(1-rhoe^2) for supply only, table 13 (1)
sigmaeta = 0.32*phic/(phic-1)*sqrt((1-rhoeta1^2*(1+rhoeta2)/(1-rhoeta2)-rhoeta2^2)); % 0.32, 0 for table 13 (1), 0.45 for table 13 (2)-(7)
deltaa = sign(phic-1)*0.02; % 0.02 for benchmark, 0 for table 13 (1), 0.05 for table 13 (2)-(7)
ecap = 1/3;

% Policy
rhor = 0.7;
rhopi = 1.5; % 1.5, 2 for Table 12 (5) and Table 13 (3)-(7)
rhoy = 0.1;
sigmam = 0.0032;
pi_ss = log(1.03);
yn1_ss = -log(delta*exp(mu)^(-1/psi))+pi_ss;
%
a1k = (exp(mu)-1+deltak)^(1/xik);
a2k = 1/(1-xik)*(exp(mu)-1+deltak);

% foreign sdf
rhofx = 1-0.035^2/(2*0.4162^2); % assume an annual vol of fx is 7%
xifx = sqrt(0.4162^2*(1-rhofx^2));

% Computation of Steady State
l_ss = 1/3 ;
k_ss = ((delta^(-1)*exp(1/psi*mu)-(1-deltak))/(alpha*(1-1/nu)*exp((1-alpha)*a0)))^(1/(alpha-1))*l_ss;
y_ss = exp((1-alpha)*a0)*k_ss^alpha*l_ss^(1-alpha);
i_ss = k_ss*(exp(mu)-1+deltak);
cc_ss = y_ss-i_ss;
ce_ss = cc_ss;
pe_ss = (1-alphac)/alphac;
c_ss = (alphac*cc_ss^((phic-1)/phic)+(1-alphac)*(ce_ss)^((phic-1)/phic))^(phic/(phic-1));
m_ss = delta*(exp(mu))^(-1/psi);
r_ss  = -log(m_ss) ;
rf_ss  = -log(m_ss) ;
Fk_ss = (1-1/nu)*alpha*y_ss/k_ss;
Fl_ss =(1-1/nu)*(1-alpha)*y_ss/l_ss;
Jk_ss  = -(alpha/nu)/k_ss;
Jl_ss  = -((1-alpha)/nu)/l_ss;
w_ss =  Fl_ss;
ph_ss = (alphac^phic+(1-alphac)^phic*pe_ss^(1-phic))^(1/(1-phic)); // 18
omega = (nu-1)/nu*((lbar-l_ss)/c_ss)^(1/xil)*(1-alpha)*y_ss/(ph_ss*l_ss); 
Phik_ss  = i_ss/k_ss;
Phik1_ss  = 1;
qk_ss = 1;
philev = 2;
pd_ss = log(m_ss*(exp(mu))/(1-m_ss*(exp(mu))));
rhopd = exp(pd_ss)/(1+exp(pd_ss)); 

%----------------------------------------------------------------
% 5. Model
%----------------------------------------------------------------
model; 
% Household (6)
uc^(1-1/psi) = 1-delta+delta*(Euc)^(1-1/psi); // 1, utility
Euc^(1-gamma) = (uc(+1)*exp(dcl(+1)))^(1-gamma); // 2, continuation value scaled by composite consumption
mh   = delta*exp(-1/xil*dc)*exp(dcl*(1/xil-1/psi))*(uc*exp(dcl)/Euc(-1))^(1/psi-gamma); // 3, SDF in consumption basket
c^((phic-1)/phic) = alphac*(cc)^((phic-1)/phic)+(1-alphac)*exp((phic-1)/phic*eta)*ce^((phic-1)/phic); // 4, consumption basket
m = mh * ph(-1)/ph; // 5, SDF in core goods
w/ph = omega*((lbar-l)/c)^(-1/xil); // 6, labor supply

% Production and the market clearing (2)
y =  k(-1)^(alpha) * l^(1-alpha) * exp((1-alpha)*a); // 7, core production function, k_t is k(-1)
cc + i +  Phip =  y; // 8, market clearing of core goods

% Physical capital (8)
k*exp(dn) = (1-deltak)* k(-1)+ Phik* k(-1); // 9, physical capital law of motion
Phik = a1k/(1-1/xik)* (i/k(-1))^(1-1/xik)+a2k; // 10, adjustment cost of physical capital
exp(rk) =  1/qk(-1)*( Fk-lambda*Jk+ qk*(1-deltak- Phik1*i/k(-1) + Phik)); // 11, Return to physical capital
Fk = alpha*(1-1/nu)* y/k(-1); // 12, auxiliary, Fk
Jk = -(alpha/nu) / k(-1); // 13, auxiliary, Jk 
Phik1 = a1k* (i/k(-1))^(-1/xik); // 14, auxiliary, Phik'
m(+1)*exp(rk(+1))  =1; // 15, Euler equation for physical capital return
qk*Phik1 = 1;  // 16, Q-theory for physical capital
  
% Labor (3)
Fl- w - lambda*Jl = 0; // 17, FOC with labor
Fl = (1-alpha)*(1-1/nu)*y/l; // 18, auxiliary, Fl
Jl = -((1-alpha)/nu)/ l; // 19, auxiliary, Jl

% Price setting (4) 
-Phip1 - m(+1)*Phip2(+1)*exp(dn) + lambda = 0; // 20, price setting, pi
Phip = phip/2*( exp(pi-pi_ss)-1)^2*y; // 21, auxiliary, Phip
Phip1 = phip*(exp(pi-pi_ss)-1)*exp(pi-pi_ss)*y; // 22, auxiliary, Phip1
Phip2 = -phip*(exp(pi-pi_ss)-1)*exp(pi-pi_ss)*y; // 23, auxiliary, Phip2

% Monetary Policy (3)
m(+1)/exp(pi(+1))*exp(yn1)= 1; // 24, Euler equation for nominal short rate
yn1  - yn1_ss = rhor*(yn1(-1) - yn1_ss) + (1-rhor)*rhopi*(pi-pi_ss) + (1-rhor)*rhoy*log(y/y_ss) + sigmam*em; // 25, Taylor rule, short rate
m(+1)*exp(pihnum(+1)-pi(+1))*exp(rf) = 1; // 26, real interest rate

% Prices
pe = (1-alphac)/alphac*exp((1-1/phic)*eta)*(cc/ce)^(1/phic); // 27, energy price
ph*c = cc + pe*ce; // 28, headline price
pie = log(pe) - log(pe(-1)) + pi; // 29, energy inflation
pih = log(ph) - log(ph(-1)) + pi; // 30, headline inflation
exp(pihnum) = exp(pi)*(1+pe)/(1+pe(-1)); // 31, numeraire inflation
epe = log(pe(+1)) - log(pe); // 32, expected energy price change (in core)

% Exogenous processes
a = (1-rhoa)*a0 + rhoa*a(-1) + phia*sigmaa0*exp(sigmaa(-1))*en ; // 33, Exogenous TFP process
dn = (1-rhon)*mu + rhon*dn(-1) + phin*sigmaa0*exp(sigmaa(-1))*en; // 34, permanent shock
log(ce) = (1-rhoe)*log(ce_ss) + rhoe*log(ce(-1)) + sigmae*sqrt(1-rhoe^2) * ee; // 35, energy supply shock
etalag = eta(-1); // 36, lag energy demand (only useful when AR(2))
eta = rhoeta1*eta(-1) + rhoeta2*etalag(-1) + sigmaeta*eeta + deltaa*en; // 37, energy demand shock
sigmaa = rhosigma*sigmaa(-1) + sigmasigma*esigma - rhosigma_a*sigmaa0*en; // 38, volatility


% asset price (21)
d = y-w*l-i-Phip; // 39, dividend in core
dagg = d + ecap*pe*ce; // 40, aggregate dividend in core

dd = log(d) - log(d(-1)) + dn(-1); // 41, real dividend growth in core
ddr = dd - pihnum + pi; // 42, dividend growth in numeraire
ddn = ddr + pihnum;  // 43, nominal dividend growth

ddagg = log(dagg) - log(dagg(-1)) + dn(-1); // 44, real aggregate dividend growth in core
ddaggr = ddagg - pihnum + pi; // 45, aggregate dividend growth in numeraire
ddaggn = ddaggr + pihnum; // 46, nominal aggregate dividend growth

exp(pd) = m(+1)*(exp(pd(+1))+1)*exp(dd(+1)); // 47, Euler equation in core price-dividend ratio (pd ratio is unit free)
exp(rd) = (exp(pd)+1)*exp(dd)/exp(pd(-1)); // 48, core equity return expressed in price-dividend ratio (in core)
exp(rdr) = (exp(pd)+1)*exp(ddr)/exp(pd(-1)); // 49, core equity return expressed in price-dividend ratio (in numeraire)

exp(pdagg) = m(+1)*(exp(pdagg(+1))+1)*exp(ddagg(+1)); // 50, Euler equation for aggregate stock market
exp(rdagg) = (exp(pdagg)+1)*exp(ddagg)/exp(pdagg(-1)); // 51, aggregate equity return (in core)
exp(rdaggr) = (exp(pdagg)+1)*exp(ddaggr)/exp(pdagg(-1)); // 52, aggregate equity return (in numeraire)
exp(rdaggn) = (exp(pdagg)+1)*exp(ddaggn)/exp(pdagg(-1)); // 53, aggregate nominal equity return

rm = log(exp(rf(-1))+philev*(exp(rdr)-exp(rf(-1)))) + sigmad*ed; // 54, market equity return (levered)
rmagg = log(exp(rf(-1))+philev*(exp(rdaggr)-exp(rf(-1)))) + sigmad*ed; // 55, aggregate equity return (levered)
rmex = rm(+1)-rf; // 56, expected excess return of core equity
rmexagg = rmagg(+1)-rf; // 57, expected excess return of the aggregate equity return

Em = log(m(+1)); // 58
m(+1)/exp(pi(+1))*fen = m(+1)*pe(+1); // 59

% The next few equations define bond yields from 1-quarter to 20-quarter (19)
exp(-2*yn2) = m(+1)/exp(pi(+1))*exp(-yn1(+1)) ; // 60
exp(-3*yn3) = m(+1)/exp(pi(+1))*exp(-2*yn2(+1)) ; // 61
exp(-4*yn4) = m(+1)/exp(pi(+1))*exp(-3*yn3(+1)) ; // 62
exp(-5*yn5) = m(+1)/exp(pi(+1))*exp(-4*yn4(+1)) ; // 63
exp(-6*yn6) = m(+1)/exp(pi(+1))*exp(-5*yn5(+1)) ; // 64
exp(-7*yn7) = m(+1)/exp(pi(+1))*exp(-6*yn6(+1)) ; // 65
exp(-8*yn8) = m(+1)/exp(pi(+1))*exp(-7*yn7(+1)) ; // 66
exp(-9*yn9) = m(+1)/exp(pi(+1))*exp(-8*yn8(+1)) ; // 67 
exp(-10*yn10) = m(+1)/exp(pi(+1))*exp(-9*yn9(+1)) ; // 68
exp(-11*yn11) = m(+1)/exp(pi(+1))*exp(-10*yn10(+1)) ; // 69
exp(-12*yn12) = m(+1)/exp(pi(+1))*exp(-11*yn11(+1)) ; // 70
exp(-13*yn13) = m(+1)/exp(pi(+1))*exp(-12*yn12(+1)) ; // 71
exp(-14*yn14) = m(+1)/exp(pi(+1))*exp(-13*yn13(+1)) ; // 72
exp(-15*yn15) = m(+1)/exp(pi(+1))*exp(-14*yn14(+1)) ; // 73
exp(-16*yn16) = m(+1)/exp(pi(+1))*exp(-15*yn15(+1)) ; // 74
exp(-17*yn17) = m(+1)/exp(pi(+1))*exp(-16*yn16(+1)) ; // 75
exp(-18*yn18) = m(+1)/exp(pi(+1))*exp(-17*yn17(+1)) ; // 76
exp(-19*yn19) = m(+1)/exp(pi(+1))*exp(-18*yn18(+1)) ; // 77
exp(-20*yn20) = m(+1)/exp(pi(+1))*exp(-19*yn19(+1)) ; // 78
slope = yn20 - yn1; // 79

% Growth rate
dc = log(c) - log(c(-1)) + dn(-1); // 80, consumption growth, notice that dn = log(N_t+1)-log(N_t)??
dcc = log(cc) - log(cc(-1)) + dn(-1); // 81, core consumption growth
dce = log(ce) - log(ce(-1)) + dn(-1); // 82, energy consumption growth
dz = a - a(-1) + dn(-1); // 83, TFP growth
Edz = dz(+1); // 84, expected TFP growth
cl = (1/(1+omega)*c^((xil-1)/xil)+omega/(1+omega)*(lbar-l)^((xil-1)/xil))^(xil/(xil-1)); // 85, bundle of consumption and leisure
dcl = log(cl) - log(cl(-1)) + dn(-1); // 86, composite consumption growth
dcr = dc + pih - pihnum; // 87, consumption growth in numeraire
dy = log(y+ce*pe) - log(y(-1)+ce(-1)*pe(-1)) + dn(-1) - pihnum + pi; // 88, output growth in numeraire
di = log(i) - log(i(-1)) + dn(-1) - pihnum + pi; // 89, investment growth in numeraire

% Cash flow and discount rate decomposition
cf = ddaggr(+1) + rhopd*cf(+1); // 90, cash flow component of pd ratio
dr = rdagg(+1) + rhopd*dr(+1); // 91, discount rate component of pd ratio
cfn = ddaggn(+1) + rhopd*cfn(+1); // 92, cash flow component (nominal) of pd ratio
drn = rdaggn(+1) + rhopd*drn(+1); // 93, discount rate component (nominal) of pd ratio
drrf = rf  + rhopd*drrf(+1); // 94, risk-free component of DR rate
dry = yn1 + rhopd*dry(+1); // 95, risk-free (nominal)component of DR rate
Erd = rdaggr(+1)  ; // 96, expected return (real, in numeraire)
Erdn = rdaggn(+1); // 97, expected nominal return


% foreign SDF
% log(mf) = rhofx*(log(m/exp(pi)))+xifx*efx;
% mf(+1)*exp(rf_for) = 1; 
end;
%----------------------------------------------------------------
% 6. Steady State
%----------------------------------------------------------------
steady_state_model ; 
uc = ((1-delta)/(1-delta*exp(mu)^(1-1/psi)))^(1/(1-1/psi)); // 1
Euc = uc*exp(mu) ; // 2
m = delta*(exp(mu))^(-1/psi); // 3
mh = m; // 4
y = y_ss; // 5
cc = cc_ss; // 6
ce = ce_ss; // 7
c = c_ss; // 8
i = i_ss; // 9
k = k_ss; // 10
dn  = mu; // 11
l = l_ss; // 12
rk = rf_ss; // 13
pi = pi_ss; // 14
pih = pi_ss; // 15
pie = pi_ss; // 16
pihnum = pi_ss;
pe = pe_ss; // 17
ph = ph_ss; // 18
w = w_ss; // 19
yn1  = yn1_ss ; // 20
lambda  = 0; // 21
qk  = 1; // 22
a = a0; // 23
eta = 0; // 24
dc  = mu; // 25
dcl  = mu; // 26
dcr = mu;
dcc = mu; // 27
dce = mu; // 28
dz  = mu; // 29
Edz = mu; // 30
Phik  = Phik_ss; // 31
Phik1  = 1; // 32
Fk  = Fk_ss; // 33
Fl  = Fl_ss; // 34
Jk  = Jk_ss; // 35
Jl  = Jl_ss; // 36
Phip  = 0; // 37
Phip1 = 0; // 38
Phip2 = 0; // 39
rf = rf_ss; // 40
cl = (1/(1+omega)*c_ss^((xil-1)/xil)+omega/(1+omega)*(lbar-l_ss)^((xil-1)/xil))^(xil/(xil-1)); //41
d = y - w*l - i; // 42
dd = mu; // 43
ddn = mu + pi;
ddr = mu;
pd = pd_ss; // 44
rd = log((exp(pd)+1)*exp(mu)/exp(pd)); // 45
dagg = d + ecap*pe*ce; // 46
pdagg = pd; // 47
rdagg = rd; // 48
rdr = rd;
rdaggr = rd;
ddagg = mu; // 49
ddaggn = ddagg+pi; // 50
ddaggr = ddagg;
rdaggn = rdagg+pi; // 51

yn2 = yn1; // 52
yn3 = yn1; // 53
yn4 = yn1; // 54
yn5 = yn1; // 55
yn6 = yn1; // 56
yn7 = yn1; // 57
yn8 = yn1; // 58
yn9 = yn1; // 59
yn10 = yn1; // 60
yn11 = yn1; // 61
yn12 = yn1; // 62
yn13 = yn1; // 63
yn14 = yn1; // 64
yn15 = yn1; // 65
yn16 = yn1; // 66
yn17 = yn1; // 67
yn18 = yn1; // 68
yn19 = yn1; // 69
yn20 = yn1; // 70
slope = 0; // 71
dy = mu; // 72
di = mu; // 73
rm = log(exp(rf)+philev*(exp(rd)-exp(rf))); // 74
rmagg = log(exp(rf)+philev*(exp(rdagg)-exp(rf))); // 75

rmex = rm - rf; // 76
rmexagg = rmagg - rf; // 77
Em = log(m); // 78
fen =  pe*exp(pi_ss); // 79
cf = mu/(1-rhopd); // 80
dr = r_ss/(1-rhopd); // 81
cfn = (mu+pi)/(1-rhopd);
drn = (r_ss+pi)/(1-rhopd);
drrf = r_ss/(1-rhopd); // 82
dry = yn1/(1-rhopd); // 83
Erd = r_ss; // 84
Erdn = r_ss+pi; // 85
sigmaa = 0; //86
epe = 0; // 87
etalag = 0; // 88
end;
steady;

shocks;
  var en = 1;
  var em = 1;  
  var ee = 1;
  var eeta = 1;
  var ed = 1;
  var esigma = 1;
end;
 
stoch_simul(irf = 0, order = 3 ,periods = 20000, pruning);



